sample.a.mean<-function(n,mu,sigma){return(rnorm(1,mu,sigma/sqrt(n)))}

plot.some.means<-function(n,sigma,mus,number.of.means,height,low,high,add=FALSE){
  J <- length(mus)
  if(!add)plot(NULL,NULL,xlim=c(low,high),ylim=c(-1,1),yaxt="n",
               xlab="sample.mean",ylab="")
  abline(h=height)
  for(i in 1:number.of.means){
    for(j in 1:J){
      cols=c("red","blue","green")
      x <- sample.a.mean(n,mus[j],sigma)
      Sys.sleep(1)
      points(x,height,col=cols[j],cex=2,lwd=2)
    }
  }
}

plot.some.means(25,15,c(100,100,100),10,0,80,120)
plot.some.means(25,15,c(90,100,110),10,0.5,80,120,add=TRUE)